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Abstract 

We show how to increase the order of one-dimensional discrete 
gradient numerical integrator without losing its advantages, such as 
exceptional stability, exact conservation of the energy integral and 
exact preservation of the trajectories in the phase space. The accuracy 
of our integrators is higher by several orders of magnitude as compared 
with the standard discrete gradient scheme (modified midpoint rule) 
and, what is more, our schemes have very high accuracy even for large 
time steps. 
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1 Introduction 



In this paper we introduce and develop discrete gradient schemes of high 
order. Discrete gradient schemes are useful tools for numerical integration 
of many-body dynamical systems [JJ El El El E] . They preserve exactly (up 
to round-off errors) both the total energy and angular momentum. More 
recently discrete gradient methods have been extended and developed in the 
context of geometric numerical integration j6]. Quispel and his coworkers 
constructed numerical integrators preserving all integrals of motion of a given 
system of ordinary differential equations [7J El El [10]. Similar ideas were 
applied in molecular dynamics simulations of spin liquids [IT] . 

In general, geometric numerical integrators are very good in preserving 
qualitative features of simulated differential equations but it is not easy to 
enhance their accuracy (in the same time preserving their geometric prop- 
erties). Symplectic algorithms can be improved using appropriate splitting 
methods [121 [T31 [HI [T51 [TBI El- Our research is concentrated on improv- 
ing the efficiency of the discrete gradient method (which is not symplectic) 
without loosing its outstanding qualitative advantages. Results reported in 
earlier papers are very promising [HI [19] . 

In this paper we present a further essential improvement of our approach 
constructing discrete gradient schemes of any prescribed order iV for one- 
dimensional Hamiltonian systems of the form: 



where V(x) is a potential, and the dot and the prime denote differentiation 
with respect to t and x, respectively. In this case the discrete gradient method 
reduce to the so called modified midpoint rule: 




(1) 



h 



n 




(2) 



Pn+l - Pn 
h 



n 



where h is the time step. 



2 



2 Discrete gradient schemes of iVth order 

We consider the following family of nonstandard numerical schemes (param- 
eterized by a single function 5): 

— 2 — = g ^ n+1 + • 

(3) 

Pn+l ~ Pn _ V(x n+ i) ~ V(x n ) 
& Xn+l X n 

where 5 can depend on any variables and parameters, including h,x n ,p n , 
x n+ i,p n+ \. One can easily prove that the total energy is preserved, i.e., 

-p 2 n + V{x n ) — E — const , (4) 

for any choice of the function 5. This is an essential generalization of the 
well known case 5 = h. 

In our recent papers [TBI HI 120] we consider 5 of the form 

( J=l t aa^ > u = y/V^, (5) 

where x may depend on x n , x n+1 but usually does not depend on h. Taking 
X — Xq } where V'(xq) = 0, we get the modified discrete gradient scheme 
(MOD-GR) [18J. Then, x = x n and x = \{x n + x n+ i) yield locally exact 
discrete gradient scheme (GR-LEX) and its symmetric modification (GR- 
SLEX), see [19j [20] . These three numerical methods are of second, third, 
and fourth order, respectively [T9] . 

In the present paper we will show that the family of numerical integra- 
tors of the form ([3]) contains numerical schemes of any order. The explicit 
formulae will be presented up to the order 11. 

The system (jHJ) (where x n = x and p n = p are given and 5 n = 5 is a 
small parameter) implicitly defines x n+ \ and p n +i- Therefore, using implicit 
differentiation, we can write down the corresponding Taylor series: 

x n+1 = x + P 5- \V5 2 - \ P V"5 3 + ^ (3V'V" - 2p 2 V") <5 4 + 0{5 5 ) , 

p n+1 =p -v>5- \pV"5 2 + i (WV" - 2V"'p 2 ) 5 3 (6) 

-i (4#T"' + 3p(V") 2 - p 3 V^) 5 4 + 0{5 5 ) . 
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Now, we assume that x n+ \ and p n+ \ coincide with the exact solution up to 
the order N, i.e., their expansion in Taylor series have at least N first terms 
identical with the Taylor series (|2T|) . Then, we compute the first N terms of 
the Taylor series of S using the first equation of Q, i.e., 

^ 2(x n _)_i x n ) . . 

Pn+1 +Pn ' ' 

The resulting polynomial of iVtli order is denoted by Sjy and its coefficients 
are denoted by a k , i.e., 

N N 

5 N = S N (x,p, h) = 2J a k (x,p)h k — h + 2_. dk(x,p)h k (8) 

k=l k=3 

where a\ = h, a 2 = 0, and a k (for k > 3) are polynomials with respect to 
p with coefficients depending on x through derivatives of V. Denoting by a 
subscript the differentiation with respect to x (and using abbreviations like 
V^x = V xxxx ) we present some number of coefficients a k in an explicit form: 

0>3 = JxVxx , «4 = ^pVm , (9) 

a 5 = ^{2V x 2 x -4V x V xxx + 3p 2 V 4x ) , (10) 
X 

a 6 = — ((5V XX V XXX - 15V x V 4x )p + W 5xP 3 ) , (11) 



<i 7 



(a 70 + a 72 p 2 + a 1A p^) , (12) 



20160 



«8 = (asip + a S3 p 3 + a 85 p 5 ) , 

725760 

aio = » ot -L nn ( a ioiP + cti03P 3 + ai05P 5 + a W7 p 7 ) , 
7257600 

an = 1596 g 72 Qo ( ano + au2p2 + ai14 ^ 4 + anep6 + a n8P 8 ) > 



(13) 
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where the coefficents a jfc , a^ m depend on x through derivatives of V, namely: 
a 70 = 17V x l + 45V x 2 V Ax - UV X V XX V XXX , 

a 7 2 = 20V? XX - 12V XX V AX - 72V x V 5x , (14) 
a 74 = 10V 6x , 

08i = 210^ - A2V X V X 2 XX + 63V x 2 V 5x , 

a 83 = UV xxx V 4x - 21V xx V 5x - 35V x V 6x , (15) 
a 85 = 3V 7x , 

a 90 = 62\£ - 228V X V 2 X V XXX + IQ8(V X 2 V X 2 XX - V*V 5x ) 
+ 90V x 2 V xx V ix , 



a 9 2 = 75V XX V 2 XX + 8lV x 2 x V Ax - 462V x V xxx V 4x 

+ 360V x V xx V 5x + A20V 2 V 6x , 
a 94 = 42VI - 120(^^6, + V x V 7x ) , 

a 96 — 7Vg x , 

aioi = 460V2 X V XXX - U70V X V XX V2 XX - Q30V x V 2 x V Ax 
+ 2385V;V^y 4 x - 9^V 2 V xx V 5x - 12Q0V^V 6x , 

a 103 = 150V X 3 XX + 15^^^ - 94514 V£ - ^6V x V xxx V 5x 
+ 4831^5* + 17851^1^ + 1080^y 7s , 

a 105 = 126^5* - lUV 3x V 6x - 26lV xx V 7x - 189V x V 8x , 

a 107 = 8Vg x , 



(16) 



(17) 
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a uo = 138214 - 64481414^3* + 414014^ 14* 

+ mm ;'i u,- + 8^oi ;;\ ,,1:;,. + 7368i/ 2 i4^ 

+ 3150^^ , 
ana = 3240t£t& - 648014K?* + 696^14* 

+ 48721414*^14, + 15660K 2 V£ - 914414^5* 

+ 11988^14^ - 21000^14^ - iosoo^v^ , 

(18) 

On4 = mOV&V^ - 18():U; r ,\r. - 86761414* Vs* 
+ 1260^, V 3x V 5x + 4770V^F ta + 306014^6* 
+ 114001414*14* + 47251/V 8x . , 

aii 6 = 336K?. - 780V 3x V 7x - 980^^ - 56014^9* 

+ 12014^ , 

Qii8 = 1814o* • 

The numerical scheme ([3]), where 5 = 5n is defined by (jSJ), will be denoted 
by GR-Af. The method GR-iV is of (at least) iVth order. We point out that 
a\ — 1, a,2 = 0. It means that the methods GR-1 and GR-2 coincide with the 
discrete gradient method (GR) given by (jSJ) (in particular, GR-1 is of 2nd 
order). Actually, if the potential V is linear in x, then any method GR-N is 
exact (i.e., its order becomes infinite). 

3 Explicit Taylor schemes of Nth order 

In this section we derive explicit numerical schemes of any order, using stan- 
dard Taylor expansions. We expand x(t + h) and p(t + h) in Taylor series: 

k=0 k=0 
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where all derivatives can be replace! by functions of x,p using (P) and its 
differential consequences (e.g., p = — V"(x)x = —V"(x)p). Thus we get 

x(t + h)=x+ph- \V'h 2 - \pV"h 3 + ± (V'V" - V"'p 2 ) h 4 + 0{h 5 ) 

p(t + h)=p- V'h - \pV"h 2 + \ {V'V" - V'"p 2 ) /i 3 (20) 

+± (3pV'V"' + p{V") 2 - p 3 V^) h 4 + 0{h 5 ) . 

Therefore, the Taylor expansion can be represented in the form 

00 h k 00 h k 

x (t + h) = ^— b k (x,p) , p(t + h) = ^— c k (x,p) , (21) 

fc=0 ' fc=0 

where b k = 4^x, c k = JxP and we compute these derivative using (pQ). For 
instance, b = x, bi = x = p and 62 = i = P = —V'(x). In general, 

d , dbk . dbk . dbk Trl/ x db k , . 

at ox op ox op 

Then, p = x implies 

Cfc = ^-bk = h+i ■ (23) 
at 

The coefficients b k (k = 1, 2, ... , 11), computed recursively from (122]) . read 
b = x , h=p , b 2 = -V x , 63 = -pKx , 

^4 = V X V XX p V^aja, , 

h = p(V 2 x + 3V X V XXX ) - p 3 V 4x , 



b 6 = -W 2 V XXX - V X V 2 X + p 2 (5V xx V xxx + 6V x V 4x ) - p 4 V 5 



(24) 



5x 



h = -P(V X 3 X + 18V X V XX V XXX + 15V 2 V 4x ) 
+ p 3 (^v 2 xx + nv xx v 4x + 10V x V 5x ) - p 5 V ex , 

^8 = V X V X 3 X + M;V,,V,,, + 15V x 3 V 4x 

- p 2 (21V x 2 x V xxx + 33V X V* XX + 81V x V xx V 4x + i5V^V 5x ) (25) 
+ p 4 (21V 3x V 4x + 21V xx V 5x + 15V X V 6X ) - p e V 7x , 
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b 9 = p(V x 4 x + 81V x V 2 x V 3x + 84V 2 V 2 X + 225V 2 V xx V 4x + 105V x 3 V 5x ] 
- P 3 (75V xx V 3 l + l02V 2 x V 4x + 231V x V 3x V 4x + 255V x V xx V 5x 
+ 105^1^) + p\2lV 2 x + 4214.14, + 3614*14* + 21V x V 7x ) 
-P 7 V 8X , 

ho = -{V.VL + 81V 2 V x 2 x V 3x + 8AV X % 2 X + 225V*V xx V 4x 

+ 105V x 4 V 5x ) +p 2 (85VlV 3x + 555V X V XX V 2 X + 837V x V 2 x V 4x 

+ 1086KV 3 ^4x + 1305^14*^5* + 420VjV aE ) 

- p A (75Vl + 585V xx V 3x V 4x + 33QV X V 2 X + 357V 2 x V bx 

+ 6961^* 1<5* + M5V x V xx V 6x + 210V 2 V 7x ) + p 6 (8iV 4x V 5x 

+ 7814,14* + 5714*14* + 281414.) - P S V 9X , 

fen = -PiVL + 336V x V x 3 x V 3x + 152AV 2 V XX V 3 2 X + 2AWV 2 V 2 x V 4x 
+ 256514 3 y 3x y 4 x + 325514 3 14 a: 14 c + 94514V 6:c ) 
+ p\8\W 2 x V 2 x + 8551414. + 92214L14, + 2AWV 2 V 2 X 
+ \87r>\[,.\;,\\., + 5175V 2 V 3x V 5x + 5U5V 2 V xx V 6x 
+ 729QV x V xx V 3x V 4x + 126014 3 y 7a: ) - p 5 (810V 2 x V 4x 

+ 92iv xx v 2 x + i9!).-,i;,,\:,,\; v + 18721414.14. + 100214414 

+ 18091414*14* + 14071414*14* + 378V 2 V 8x ) + p 7 (8AV 2 



6x 

5* 



+ 16214,14. + 13514,14. + 8514*14* + 361414,) - p 9 V, 



10. 



(26) 



(27) 



(2? 



Thus for any fixed iV we obtained the following explicit numerical scheme, 
denoted TAY-iV (the Taylor scheme of iVth order), 



x n+l = /J T[ h(x n , Pn) , Pn+1 = Tf C k{x n ,Pn) , (29) 

k=0 ' k=0 



where bk and Ck are defined by ( 1221) . ( 123]) and, in particular cases, by ( |24l) . 
( I25i) . ( 1261) . f l27l) and ( |28|) . Explicit integrators Tay-A^ will be used for compar- 
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ison with discrete gradient methods of high order. Moreover, they are good 
candidates for predictors when gradient methods ([3D are used as correctors. 

4 Numerical experiments 

In our recent papers we compared several discretizations of the simple pen- 
dulum equation (V(x) = — cosx) with a special stress on the long-time 
behaviour, see jTH [19]. Locally exact discrete gradient schemes (GR-LEX 
and GR-SLEX) turned out to be the best. In some tests their accuracy was 
better by several orders of magnitude in comparison to standard methods 
like leap-frog, implicit midpoint rule or the discrete gradient method (GR). 
GR-LEX and GR-SLEX yield rather similar results and in this section we 
confine ourselves to GR-LEX only. 

We are going to compare GR-LEX with algorithms of higher order in- 
troduced in the present paper, i.e., GR-iV and TAY-iV. The accuracy of 
these schemes was tested mainly for the simple pendulum, but other poten- 
tials yield similar results. We present some data for the Morse potential, 
V(x) = \e~ 2x — e.~ x , see Fig. [2j In both cases the exact solution is known. 
For simplicity we always assume the initial position at the stable equilib- 
rium, i.e., xq = 0. The details of numerical computations of the period are 
explained in [TS] and iteration procedures are described and discussed in [Tj5] 
(we apply the fixed point method and the Newton method, and iterate until 
the acuracy 10~ 16 is obtained). We point out that Sn given by (jSj) depends 
on x n ,p n and does not depend on x n+ i,p n+ i). It means that 5n is evaluated 
only once at every step. 

4.1 Global error 

Fig. [1] and Fig. [2] show the dependence of the global error of the numerical 
solutions on the time step (the global error was evaluated at t — 1207^). 
GR-3 yields almost the same results as GR-LEX. They are better than GR 
by several orders of magnitude. GR-iV (for iV ^ 5) are more accurate than 
GR-LEX by several orders of magnitude. We point out that the schemes 
GR-iV are very accurate for large time step. Actually, for small time steps 
(say, h < 0.1) the accuracy of GR-7 and GR-11 almost does not depend on 
h (actually, it even slightly decreases for smaller h). TAY-10 becomes less 
accurate than GR-7 and TAY-5 for larger h. 
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Theoretically all gradient schemes preserve exactly the energy but, of 
course, round-off errors cause some small inaccuracy, see Fig. [3J We see, that 
the energy error accumulates slowly, almost linearly but with a very very 
small slope: for t w 300 000 we have AE w 10~ 12 . 

4.2 Stability and relative error of the period 

All gradient schemes have extremaly stable period of oscillations. The sta- 
bility of the discrete gradient scheme (GR) was tested in detail in [18]. Other 
gradient schemes follow the same pattern. Fig. H] compares the average pe- 
riod (more precisely: T avg (N, 20), see [IB]) of numerical solutions produced 
by GR-7 and TAY-10. If t is not very large, then in both cases the average 
period oscillates around the exact value T t h- For longer times we clearly see 
that TAY-10 becomes less and less exact, see Fig. HI while GR-7 oscillates 
exactly in the same way, even for very very long times, e.g., t ~ 30 000 000 
at Fig. EJ 

Fig. [H] and Fig. [7J illustrate the relative error of the period. More precisely, 
we consider T avg (0, 100, 200) (similarly as in |19j), for details see [18], p. 11 
(roughly saying, we consider the first 200 periods making some averaging). 
Then, we compare the results with the exact period T t h- 

Fig. [6] presents the dependence of the relative period on the time step h. 
We see that GR-7 yields excellent results (better by 3-4 orders of magintude 
than GR-LEX). The accuracy of TAY-10 and GR-11 is (for p = 1.95 and 
h < 0.3) almost the same. The accuracy of Taylor schemes becomes relatively 
lower for greater h. 

Increasing the order of GR-iV for small h we increase the accuracy but 
only to some extent, see Fig. [7] Indeed, for h = 0.02 schemes GR-11 and 
TAY-10 yield practicaly the same accuracy as GR-7, i.e., 10~ 13 for oscillating 
motions and 10~ 9 for rotating motions, with exception of the region po ph 2, 
where the accuracy is lower for any numerical scheme. For po < 2 (oscil- 
lations) GR-7 is more accurate than GR by 7-9 orders of magnitude. For 
small po also GR-LEX and TAY-5 attain such high accuracy. For p > 2 the 
scheme GR-LEX produces almost the same results as GR-3, and both are 
less accurate than TAY-5 (GR-LEX is more and more accurate for decreas- 
ing p Q ). We point out that gradient schemes produce very stable results (i.e., 
the picture presented at Fig. [JJis time-independent. The accuracy of Taylor 
schemes decreases with time, see Fig. HJ 
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4.3 Neighbourhood of the separatrix 

The neighbourhood of the separatrix (po ~ 2 for the simple pendulum) is 
most difficult to be simulated numerically. The discrete gradient method 
(GR) turns out to be relatively good in this region, see [18], and the locally 
exact methods (GR-LEX, GR-SLEX) work almost perfectly [T9]. Here, we 
take for comparison also GR-3, GR-7, GR-11 and TAY-5, TAY-10. The Tay- 
lor schemes are much worse in this region: for h = 0.9 even TAY-10 is not 
able to reproduce the correct qualitative behaviour, see Fig. [HJ Through- 
out the first period the scheme GR yields good qualitative behaviour and is 
better than TAY-10 with a halved time step, see Fig. |HJ In the first period 
GR-3, GR-7 and GR-LEX (and also GR-SLEX and GR-iV for N > 3) pro- 
duce similar results. We point out that the exact trajectory is very close to 
the separatrix (|po — 2| = 10~ 10 ) and h is very large but, nevertheless, all 
improved discrete gradient methods simulate very accurately the motion of 
the pendulum. 

Fig. IH] shows the same situation for much longer times (t > 100 000). 
Note that the time step for TAY-10 is much smaller (h = 0.09) than the 
time step for all gradient schemes, which is very large (h = 0.9). In spite of 
that essential handicap, TAY-10 is only slightly better than GR-7 and less 
accurate than GR-11. GR-7 is more accurate than GR-LEX. 

5 Conclusions 

The numerical integrators GR-iV, described in this paper, have similar advan- 
tages as GR-LEX and GR-SLEX: they preserve exactly the energy integral 
(i.e., eq. (Jl]) holds), are extremaly stable and have very good long-time be- 
haviour of numerical solutions. They can be constructed for any prescribed 
order N. 

Therefore, modifications presented in this paper essentially improve the 
discrete gradient method (at least in the one-dimensional case) keeping all 
its advantages. Schemes GR-iV (for iV ^ 7) are much more accurate than 
GR-LEX for most choices of parameters. Only in the region of small po the 
scheme GR-LEX is comparable with discrete gradient methods of high order. 

We point out that numerical schemes ([3]), like all discrete gradient meth- 
ods, are neither symplectic nor volume-preserving. Moreover, schemes GR-iV 
are not time-reversible. Therefore, the conservation of the energy integral 
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plus high order seem to be sufficient to assure oustanding qualitative and 
quantiaive properties of these methods. 

Acknowledgments. This research work has been supported by the grant No. 
N N202 238637 from the Polish Ministry of Science and Higher Education. 
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Figure 1: Global error at t — 1201^ as a function of the time step h for the 
simple pendulum, p = 1.8 (T th = 9.122 196 55). 
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Figure 2: Global error at t — 1207^ as a function of the time step h for 
Morse potential, for p = 0.8 (T th = 10.471 975 51). 
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Figure 3: Energy error as a function of time (t = hN, h = 0.25), for the 
simple pendulum, p = 1.8 (E ex = 0.62). 
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Figure 4: Average period as a function of time (N is a number of half-periods) 
for the simple pendulum, p — 1.8 (T t h = 9.122 196 55). Dark points - GR-7, 
light points - TAY-10, solid straight line - exact period. 
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Figure 5: Average period as a function of time (N is a number of half-periods) 
for the simple pendulum, p = 1.8, scheme GR-11. Solid straight line - exact 
period (T th = 9.122 196 55). 



18 



10° 




Figure 6: Relative error of the period of the simple pendulum as a function 
of h, for p = 1.95 (T th = 11.657 585 28). 
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Figure 7: Relative error of the period for the simple pendulum as a function 
of po, for h = 0.02. 
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Figure 8: x n as a function of time (t = nh), very near the separatrix (p = 
1.999 999 999 9), h = 0.09 for TAY-5, h = 0.45 for TAY-10, h = 0.9 for 
all other discretizations. The solid line corresponds to the exact solution 
(T th = 51.596 879 14). 
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Figure 9: x n as a function of time (t = nh), very near the separatrix (p = 
1.999 999 999 9), h = 0.09 for TAY-10 and h = 0.9 for all other discretiza- 
tions. The solid line corresponds to the exact solution (T^ = 51.596 879 14). 
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